﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace Nature.Toolkit
{

    public class SpeciesMassFractionDecompositionObject
    {
        decimal[] m_y0;
        double[] m_yCurrent;
        decimal[] m_dyTotal;
        readonly int Nsp;

        public SpeciesMassFractionDecompositionObject(ReadOnlyArray<double> y)
        {
            int numberOfSpecies = y.Length;
            m_y0 = new decimal[numberOfSpecies];
            m_yCurrent = new double[numberOfSpecies];
            m_dyTotal = new decimal[numberOfSpecies];

            for (int isp = 0; isp < numberOfSpecies; ++isp)
            {
                m_y0[isp] = (decimal)y[isp];
                m_dyTotal[isp] = 0m;
            }
            
            Nsp = numberOfSpecies;

            decimal ysum = 0m;
            for (int isp = 0; isp < Nsp; ++isp)
            {
                decimal yk = (decimal)y[isp];
                if (yk < 0m) { yk = 0m; }
                else if (yk > 1m) { yk = 1m; }
                ysum += yk;
                m_y0[isp] = yk;
            }
            for (int isp = 0; isp < Nsp; ++isp)
            {
                m_y0[isp] /= ysum;
            }
        }

        public ReadOnlyArray<double> Y { get { return m_yCurrent; } }

        public void Increment(double[] dy)
        {
            decimal ysum = 0m;
            for (int isp = 0; isp < Nsp; ++isp)
            {
                decimal dyk_prev = m_dyTotal[isp];
                decimal dyk_add = (decimal)dy[isp];
                decimal dyk_total = dyk_prev + dyk_add;
                decimal yk0 = m_y0[isp];
                decimal yk = yk0 + dyk_total;                
                if( yk < 0m) 
                {
                    yk = 0m;
                    dyk_total = (-yk0);
                }
                else if (yk > 1m)
                {
                    yk = 1m;
                    dyk_total = (1m - yk0);
                }
                ysum += yk;
                m_dyTotal[isp] += dyk_total;            
            }
            for (int isp = 0; isp < Nsp; ++isp)
            {
                double dyk = (double)(m_dyTotal[isp] /= ysum);
                m_yCurrent[isp] = (double)m_y0[isp] + dyk;
            }
        }
    }
}
